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We investigate numerically the quasiparticle density of states g(E) for a two-dimensional, dis- 
ordered superconductor in which both time-reversal and spin-rotation symmetry are broken. As a 
generic single-particle description of this class of systems (symmetry class D) , we use the Cho- Fisher 
version of the network model. This has three phases: a thermal insulator, a thermal metal, and a 
quantized thermal Hall conductor. In the thermal metal, we find a logarithmic divergence in g(E) 
as E — > 0, as predicted from sigma model calculations. Finite-size effects lead to superimposed os- 
cillations, as expected from random matrix theory. In the thermal insulator and quantized thermal 
Hall conductor, we find that g(E) is finite at E = 0. At the plateau transition between these phases, 
g(E) decreases toward zero as \E\ is reduced, in line with the result g(E) ~ \E\ ln(l/|i?|) derived 
from calculations for Dirac fermions with random mass. 

PACS numbers: 73.43.-f 73.20.Fz 74.78.-w 



I. INTRODUCTION 

The study of Anderson transitions - disorder-driven 
transitions in systems of noninteracting fermions - has a 
long history in condensed-matter physics. Symmetry is a 
guide to understanding both mesoscopic behavior in each 
phase and critical behavior at a transition. According 
to the symmetry classification inherited from random- 
matrix theory, in a conventional situation three Wigner- 
Dyson symmetry classes - orthogonal, unitary and sym- 
plectic - can be distinguished, depending on whether or 
not the Hamiltonian is invariant under time reversal and 
spin rotation (see Refs. [11131100 for reviews). Above a 
lower critical dimension, systems in each symmetry class 
can have a transition between insulating and metallic 
phases. In addition, in two-dimensional systems without 
symmetry under time-reversal or parity, a second type 
of insulating phase is possible, with edge states and a 
quantized Hall conductance, and there can be a plateau 
transition between the insulator and the quantized Hall 
conductor. Long-distance properties of these systems are 
described by nonlinear sigma models. Within this frame- 
work, quantized Hall conductors appear if symmetry al- 
lows a topological term in the sigma models 

All Anderson transitions in the standard, Wigner- 
Dyson symmetry classes are expected to share certain 
general features. One is that the density of states (DoS) 
is a smooth function of energy, and hence noncritical at 
a transition. Another is universality, in the sense that, 
for a given symmetry and dimensionality, only one type 
of behavior is expected. 

It is now widely appreciated that the range of pos- 
sibilities is not exhausted by the Wigner-Dyson classes, 
and there exist additional symmetry classes. These ad- 
ditional classes are distinguished from the standard ones 



by possessing a particle-hole symmetry that selects as 
special one particular energy E in the spectrum (we 
take this to be E = 0). Moreover, their density of 
states may be singular at E = and may have dis- 
tinct behavior in each phase and at a critical point. 
While individual examples of systems from the additional 
symmetry classes have been known for many years, as 
tight-binding models with two sublattice structure;^^ 
within random-matrix theory ; 10 : 11 and as models for 
superconductors,- 12 a full classification scheme was devel- 
oped only more recently by Altland and Zirnbauer i 13 : 14 
This scheme includes, in addition to the Wigner-Dyson 
classes, systems belonging to two additional types of sym- 
metry class. One set, the chiral classes, arises in two 
sublattice systems. The other set is realized in mod- 
els for noninteracting quasiparticles in disordered super- 
conductors based on Bogoliubov-de Gennes Hamiltoni- 
ans. The properties of models belonging to these ad- 
ditional symmetry classes have been a focus of atten- 
tion in connection with understanding low-energy quasi- 
particle DoS, transport, and localization properties in 
dirty superconductors, including those with unconven- 
tional pairing Among other fea- 
tures, the additional symmetry classes allow for versions 
of the quantum Hall effect (QHE) in two-dimensional dis- 
ordered superconductors with broken time-reversal sym- 
metry. The relevant quantized conductance in these 
cases is a thermal or spin conductance, since quasiparticle 
number is not conserved by the Hamiltonian. For super- 
conductors that are invariant under spin-rotations, the 
symmetry is referred to as class C, and the version of the 
quantum Hall effect is known as the spin QHE. Without 
spin rotation invariance, the symmetry is termed class D 
and one has the thermal QHE. The Hamiltonian for a 
system of the class D has the following block structure 
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in the particle-hole spaced: 
H=( * , * 



/it , A 



(1) 



which is determined by the condition H = —a x H T a x (in 
addition to the Hermiticity H = iff). Alternatively one 
can work in a different basis, defining H = g^Hg with 
g 2 = u x - In this basis, the defining condition of class D 
becomes H = —H T , so that H is purely imaginary. It is 
the systems with this symmetry that are the subject of 
the present paper. 

A striking feature of this class is that symmetry and 
dimensionality alone are insufficient to determine behav- 
ior. At the level of the nonlinear sigma model, the reason 
is believed to be that the relevant target space has two 
disconnected pieces and that, depending on the choice of 
underlying microscopic model, it may or may not be nec- 
essary to consider configurations containing domain walls 
on which the sigma model field jumps between the two 
components ! 21 ! 24 ' 25 In the following, we investigate the 
Cho-Fisher (CF) network model 2 ^ for the thermal QHE. 
This model is generic in the sense that it displays all three 
phases possible in a two-dimensional class D system: the 
metal, the insulator, and the quantized Hall conductor. 
By contrast, in a closely related model - a fermionic ver- 
sion of the ± J random bond Ising model (RBIM) - the 
metallic phase is absent^ 

The location of the metallic phase and the two localized 
phases in the phase diagram of the CF model is indicated 
in Fig. [1] together with phase boundaries of two types, 
corresponding to metal- insulator and plateau transitions, 
respectively. This phase diagram was obtained in Ref. [23| 
from transfer-matrix calculations; let us explain why its 
structure may be quite naturally anticipated. In the ab- 
sence of disorder (p = 0) the system has a gap around 
E = for all a except for sin 2 a = \ (where the DoS 
vanishes linearly at zero energy). One thus expects that 
the system remains insulating for weak disorder (small 
p) . More carefully, analyzing the edge states in the both 
limits of uncoupled plaquettes, sin 2 a = and sin 2 a = 1 , 
one observes that there are two topologically different in- 
sulating phases. From symmetry, if there is a direct tran- 
sition between these phases, it must occur at sin 2 a = i. 
Finally, with increasing p, the disorder creates a substan- 
tial DoS at E = 0, and the system may be expected to 
undergo a transition into the metallic phase, which gener- 
ically exists in two-dimensional systems with spin-orbit 
coupling. 

Two-dimensional systems from class D have attracted 
considerable attention recently. The sigma model 
analysi s 16 ' 21 starts from a short-distance description in 
terms of a diffusive metal with a smooth DoS that is 
finite at E = 0. From a renormalization group (RG) 
analysis (perturbative in g^ 1 , the inverse dimensionless 
conductance), it is found that the diffusion constant is 
unrenormalized at leading order and that the DoS has a 
logarithmic divergence near zero energy, which is given 
in terms of the diffusion constant D and the mean free 
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FIG. 1: Phase diagram of the Cho-Fisher model as obtained 
in Ref. |23| from transfer-matrix calculations. The plane is 
spanned by the parameters sin 2 (a), the interplaquette tun- 
neling probability, and p, the concentration of vortex disor- 
der: these control the short-distance values of the conductiv- 
ity components o xy and <r xx , respectively. 
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An alternative approach to the theory of these sys- 
tems, arguably tailored to describe the plateau transi- 
tion, has been developed in Ref. l2l| by starting from a 
model of Dirac fermions with random mass and treating 
this disorder perturbativcly, in the spirit of the analy- 
sis of the Ising model by Dotsenko and Dotsenkoi^l The 
disorder-free system has a transition, driven by tuning a 
uniform mass through zero, which in the CF model lies 
at p=0, sin 2 (a)=l/2. In the vicinity of the clean fixed 
point representing this transition, the disorder strength 
g\i is marginally irrelevant. This implies for the critical 
DoS a logarithmic correction term (see Appendix [A} of 
the form 
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Clearly, since this calculation is for a system with weak, 
homogeneous disorder, its relevance for behavior in the 
CF model with dilute, strong scatterers needs to be 
tested. 

In the localized phases, past work suggests several pos- 
sibilities for the behavior of the DoS near E=0. The 
simplest approach is to imagine that the sample can be 
divided into independent regions of size set by the local- 
ization length and that the contribution of each region 
to the DoS can be obtained from random-matrix theory 
for this symmetry class^ giving finite g(E) at E=0. Al- 
ternatives are suggested by the fact that the off-critical, 
disorder-free model (p—0, sin 2 (a)^0) has a gap in the 
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DoS around E = 0. As has been well-studied in one- 
dimensional system, rare disorder configurations (termed 
Griffiths strings) may fill in this gap, generating a DoS 
that varies as a positive or negative power of \E\ near 
■E=0i 17 i 18 Behavior of this kind has been found recently 
in the network model representation of the RBIM. 28 

The purpose of this article is to present numerical stud- 
ies of the behavior of the DoS in the CF model. An 
outline is as follows. In Sec. HH we describe the model 
and our numerical methods. Section IIIII is central for 
the paper and contains our main findings. We analyze 
first (Sec. lIII A"| the low-energy DoS in the metallic phase 
(at large p) and confirm the logarithmic divergence, Eq. 
©. As an additional manifestation of metallic behav- 
ior, we find random-matrix-theory (RMT) oscillations in 
the DoS superseding this logarithmic behavior at low- 
est energies. At smaller p, one enters a localized phase 
(Sec. IIIIB[) where neither a logarithmic divergence nor 
oscillations are observed. Instead, the DoS remains fi- 
nite at zero energy. Finally, in Sec. IIII CI we turn to the 
plateau transition which occurs for low p on the self-dual 
line sin 2 (a)=l/2 (see Fig.[T]). We find in a wide energy in- 
terval that the behavior of DoS is consistent with the RG 
result, Eq. ((3]). In Sec. lIVI we summarize our findings and 
discuss directions for future work. A brief description of 
an alternative microscopic model and results obtained for 
it are given in Appendix [Bj these results are very close 
to those for the CF model, supporting the idea that the 
latter is generic. 



II. MODEL 

The CF model 2 - 6 belongs to the family of network mod- 
els first proposed in Ref. [29| for the description of the 
(conventional) quantum Hall effect. The structure of the 
network in the clean limit is shown in Fig. O The model 
describes quantum dynamics of nonintcracting particles 
living on the directed links of a square network. If (as 
here) the wave function on a link has one component, the 
state of the system as a whole is represented by an N- 
component vector, where TV is the number of links. The 
time evolution is characterized by a N x N unitary ma- 
trix U, which in the absence of disorder is governed by 
a single parameter, < a < ir/2. At every node of the 
network the particle turns right or left with a probability 
amplitude ±cosa or ±sina, respectively. The ampli- 
tude signs are shown in Fig. [5] and ensure unitarity of 
the evolution operator Tj£2- Note that each plaquette of 
the network carries half a flux quantum so that particles 
pick up a phase factor ir when moving around it. 

In the absence of disorder, the parameter a is the same 
for all the nodes; then the network is a fermionic repre- 
sentation of the clean Ising model and U can be diago- 
nalized by Fourier transform. Near the critical point and 
for small wave vectors fc, the spectrum has a Dirac form 
E = (k 2 + £- 2 y/ 2 . The gap £ _1 = \a- tt/4| vanishes at 
the critical point a — a c = 7r/4 where the DoS is linear 
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FIG. 2: Network representation of the clean Ising model. 
The plaquettes form two sublattices, C and S. The symbols 
"s" and "c" denotes the amplitudes for left and right turns, 
±sina and ±cosa, governing the evolution at the network 
nodes. Adding disorder by inserting a vortex pair corresponds 
to flipping the signs of either a pair of s or a pair of c that are 
associated with a given node. 



in energy: g(E) — \E\/2n. 

One can think of the evolution U as being generated 
by a Hamiltonian. Taking U = exp(— iH), the symmetry 
condition for class D, that H is purely imaginary, implies 
that U is real. Hence, within symmetry class D, disorder 
can be introduced into the model by allowing for node- 
to-node variation of the parameter a. This can be done 
in a variety of ways. A disorder weak in the perturbative 
sense of Refs. andl2ll can be realized by drawing a for 
every node from a distribution with a width Sa that is 
small compared to the mean value. In contrast, in the 
CF model, disorder is introduced as isolated defects by 
making the change a — > —a or a — > tt — a, for a subset 
of nodes randomly distributed with a concentration p. 
This amounts to flipping signs of either both sin a or 
both cos a associated with such a node. This procedure 
can be viewed as the insertion of two additional half- 
flux lines into two plaquettes adjacent to the node and 
belonging to the same sublattice, see Fig. [2j Note that 
the vortex pair appears with equal probability on the C- 
or S-sublattice. It is this feature that distinguishes the 
CF model from the RBIM, 26 '- 1 which is obtained if all 
the additional vortices are placed on the same sublattice. 

Our numerical analysis centers on the matrix U for a 
system with a torus geometry of size L x L. Since U 
is unitary, its L 2 eigenvalues e lEj lie on the unit circle, 
defining the energies Ej. For a given realization of dis- 
order, we compute eigenstates in a vicinity of the value 
E=0, where the special features of class D reveal them- 
selves. For this purpose, we use efficient sparse matrix 
packages ! 32 ! 33 ' 34 The procedure is performed for an en- 
semble consisting typically of 10 4 or 10 5 disorder realiza- 
tions, and the DoS is obtained as an ensemble average. 
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FIG. 3: Low-energy DoS in the metallic phase. Parameters 
(upper curves): p=0.5, a=7r/4, system sizes L—128 (squares), 
and L=256 (full circles). The straight dashed line represents 
the logarithmic asymptotics. For lowest energies, the RMT 
oscillations are clearly visible; they can be collapsed on a sin- 
gle curve, as shown in Fig. [4] For comparison, the results for 
p = 0.15 and p = 0.1 are also shown. 



III. RESULTS AND DISCUSSION 

A. Thermal metal 

Before presenting our results for the metallic phase 
and comparing them to the analytical prediction Eq. |(5J), 
we briefly recall the theoretical framework within which 
Eq. @ is obtained j 16 i 21 Using the standard procedure, 
one derives an effective field theory that has the form of 
a diffusive supersymmetric nonlinear sigma model. This 
theory is valid on energy scales E < t -1 , where r _1 is 
the elastic transport scattering rate. For the lowest en- 
ergies, E < i?Th, below the Thouless energy E^t (the 
inverse time of diffusion through the system), the the- 
ory becomes effectively zero-dimensional and reproduces 
the random-matrix theory of class D. A renormalization- 
group analysis, perturbative in the running coupling con- 
stant / (which is proportional to g^ 1 ), yields 



df 
dln£ 



(4) 



where £ denotes the ultraviolet cutoff. This implies that 
the infrared behavior of the system is governed by the 
perfect-metal fixed point, / — > 0. In other words, with 
increasing system size L, the conductance of the metal in- 
creases logarithmically, g(L) oc \nL, and hence diverges 
in the limit L — > oo, so that the perturbative RG is jus- 
tified. A similar analysis yields the RG equation for the 
second coupling constant e, whose bare value is given by 
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FIG. 4: Renormalized DoS at maximal disorder p=0.5 and on 
the symmetry line sin 2 a=l/2 for different system sizes vs the 
energy measured in units of the level spacing &l- The RMT 
result, Eq. (0, is plotted as a dashed line for comparison. 
The inset shows the logarithmic dependence of 1/L 2 8l on 
the system size L, consistent with the data of Fig. [3] A fit of 
the slope yields g f — 0.152 corresponding to D — 0.33. 
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leading to Eq. §Z§ for the DoS. Comparing Eqs. (0| and 
([5]) , one sees that the logarithmic increase in the conduc- 
tance is driven by the logarithmically divergent density 
of states, while the diffusion constant D remains nonsin- 
gular to this order. 

We now turn to the results of our numerical simula- 
tions. Fig. shows the DoS calculated at the maximal 
concentration p of flux lines, p=l/2. The data exhibit a 
logarithmic increase of the DoS over almost three decades 
in E for the larger system size, L=256. We stress that the 
increase continues to be of logarithmic form even though 
the renormalized DoS at small energies becomes much 
larger than its bare (large- E) value go— 0.1. This is a 
signature of the fact that the RG flow is toward weak 
coupling, so that the one-loop result ([2]) is valid down to 
arbitrarily low energies in the thermodynamic limit. 

At the smallest energies, we observe pronounced os- 
cillations in the DoS. These are RMT oscillations due 
to finite system size and serve as another indication of 
the fact that we are dealing with a metallic phase. To 
demonstrate the RMT origin of these oscillations, we re- 
plot these parts of DoS curves, rescaling the energy to 
the mean level spacing Sl at lowest energy for the corre- 
sponding system size. Specifically, Sl is obtained by nu- 
merically solving q{Sl) = 1/L 2 Sl as suggested by Eq. (J7]) 
below. The results are shown in Fig. 3] for six different 
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system sizes. The data collapse on a single curve, which 
shows that the (renormalizcd) level spacing 
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is indeed the only relevant energy scale in the regime 
E < E^h where the RMT is applicable. As further shown 
in Fig.|H the curve obtained agrees nicely with the RMT 
prediction, 
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up to _E/<5l~1.5-2; for larger energies, the oscillations are 
strongly suppressed. This is fully consistent with the ex- 
ponential vanishing of the RMT oscillations beyond the 
Thouless energy (see, e.g., Ref. HI). With increasing sys- 
tem size, the ratio E^/Sl increases (though only log- 
arithmically), so that the RMT range includes progres- 
sively more oscillation periods. This tendency is clearly 
seen in Fig. [4j 

For comparison, we also show in Fig. [3] the low-energy 
DoS for p = 0.15 and p = 0.1; the latter point is close to 
the expected boundary of the metallic phase, see Fig. Q] 
It is seen that when the system approaches the phase 
boundary, the logarithmic increase of the DoS disappears 
and the RMT oscillations get damped. 



B. Localized phases 

Having explored the metallic phase, we now turn to 
the phases with localized states. Fixing the interpla- 
quette coupling a at a value somewhat different from 
sin 2 (a)=l/2 and decreasing the concentration p of vor- 
tex defects, we expect (see Fig. [1]) that the system un- 
dergoes a transition at a critical value p c {a) (which is of 
the order ~ 0.1) and enters the insulator or quantum Hall 
conductor. In both phases, bulk states are localized. Fig- 
ure [5] demonstrates how this transition is reflected in the 
behavior of the DoS. In this figure, we display the evolu- 
tion of the DoS at fixed tunneling probability sin 2 a~0.58 
with decreasing defect concentration. The low-energy 
singular peak, whose magnitude increases logarithmically 
with the system size (as studied in Sec. IIII A[) . is still 
clearly seen for p=0.13 (see upper inset of Fig. [5]) but 
is absent once p falls below 0.1 (lower inset of Fig. [5]). 
Since this logarithmic singularity was a key signature of 
the thermal-metal phase, its absence signals the break- 
down of the metallic behavior and, hence, the emergence 
of a localized phase. 

For the localized phases, our main finding is that the 
DoS has a nonzero value at E=0. It is interesting to 
ask how this behavior connects with the value g(0)—0 
expected from Ref. Hi] (see Eq. [3] above) at the plateau 
transition. Consider the DoS for a close to the criti- 
cal value 7r/4, so that the localization length £ is large. 
The behavior of DoS can then be understood by using 
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FIG. 5: DoS near E = for disorder values p = 0.13, 0.08, 
0.04, and 0.02 and for two system sizes L = 64, 128 at fixed 
interplaquette coupling sin 2 (a) = 0.579. The DoS diverges 
logarithmically as E — > in the metallic phase (p = 0.13) 
and remains finite in the localized phase (other values of p). 
The results for the lowest impurity concentration, p = 0.02, 
show an oscillatory feature induced by the band structure of 
the clean system, as well strong scatter in the data at the 
lowest energies, which is due to insufficient ensemble averag- 
ing. Upper inset: Low-energy peak at p = 0.13; its amplitude 
increases with L, in agreement with Sec. IIII Al Lower inset: 
Low-energy DoS at p = 0.08. No peak at E — > is detected; 
p(E — > 0) is a constant independent of L, indicating that the 
system is in the insulating phase. Statistical noise in the lower 
inset is more pronounced than in the upper one due to the 
smallness of the DoS. 



the Dirac-fermion RG presented in Appendix |A"1 Specifi- 
cally, for energies that are not too small, behavior will be 
the same as at criticality, Eq. ([3]). However, for smallest 
energies, it is the localization length £ (rather than E) 
that will terminate the RG process. In this sense, the 
role of £ is fully analogous to that of finite system size L 
at criticality. This implies (see Appendix [AJ that g(E) 
saturates at the value 
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where to is the ultraviolet cutoff length. The energy E% 
at which the saturation takes place is 
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The low-energy saturation of the DoS in the localized 
phases which we observe in our numerical simulations is 
fully consistent with these analytical predictions. 
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FIG. 6: DoS at low energy on the self-dual line sin 2 (a) = 0.5 
for disorder concentrations p = 0.2 (o, o), 0.1 (□, x), and 0.05 
(A, +), where in each case, the first symbol is for L = 128 
and the second is for L = 256. Inset: g(E)/\E\ at p=0.05 
on a log-linear scale. The logarithmic correction is clearly 
observed, in agreement with Eq. 



Before closing this subsection, we comment briefly on 
the regions of localized phases where the interplaquette 
coupling is very weak (sin 2 a close to zero or to unity) . As 
shown recently, 2 - in this situation the DoS of the RBIM 
acquires a nonuniversal power-law singularity, li^ 1 / 2-1 
with z > 1 associated with Griffiths stringsj 17 ' 18 We ex- 
pect that the same mechanism should be operative in the 
Cho-Fisher model as well. An analysis of these parts of 
the phase diagram and of the expected Griffiths singular- 
ities is, however, outside the scope of the present paper. 



Plateau transition: 
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The phase boundary between the insulator and ther- 
mal Hall conductor is the location of the plateau tran- 
sition. From the treatment of Ref. [2l| (see also Ap- 
pendix^, RG flow on this boundary is toward the clean 
Ising fixed point. The corresponding RG result for the 
critical DoS is given by Eq. ©. To test this prediction, 
we have studied the evolution of the DoS on the self-dual 
line sin 2 a = 1/2 with decreasing concentration p of dis- 
order. Results are shown in Fig. [S] For the largest two 
values, p — 0.2 and p = 0.1, the DoS exhibits a peak 
at E = whose amplitude increases with L, which is a 
hallmark of the metallic phase (Sec. IIIIB[) . For the low- 
est value, p = 0.05, this peak is not observed anymore, 
suggesting that this point belongs to the phase bound- 
ary between the two quantum Hall localized phases, see 
Fig. [TJ Indeed, plotting g(E)/\E\ as a function of log-E 
(see the inset to Fig. [5]), we find a behavior fully consis- 
tent with the logarithmic correction predicted by Eq. ^ . 
At lowest energies, an oscillatory structure is observed in 
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FIG. 7: Evolution of the DoS at low energies with increas- 
ing system size. DoS obtained by binning the 16 lowest-lying 
eigenvalues calculated for each network operator out of an 
ensemble of typically 10 4 disorder realizations. The step- 
like structures in the DoS is a remnant of the band struc- 
ture of the clean lattice and their effect gradually dimin- 
ishes with L increasing. Bold dashed lines show a fit to 
Eq. ©, g(E)/\E\ = ci]n\E /E\. Upper panel: p=0.08, 
L=16,32,64, 128,256,512, 1024, ci =0.524, £ =0.433. The 
position of the sharp minimum occurring at small sys- 
tem sizes, L=16, 32, indicates the mean level spacing in 
these samples. Lower panel: p=0.05, L=128, 256, 512, 1024, 
ci =0.1875,S =2.865. 



the DoS curve for p — 0.05 in Fig. [6l This feature is a 
finite-size effect and is governed by the few lowest-lying 
eigenstates which inherit information on their position 
in the clean system. With increasing system size, the 
energy window for these oscillations shrinks, so that the 
DoS acquires a smooth limiting form in the thermody- 
namic limit. 

We turn now to the analysis of the small-p DoS at 
asymptotically low energies. This requires investigation 
of system sizes larger than those used in Fig. [S] We show 
in Fig. [7] a log- log plot of the DoS for the disorder concen- 
trations p — 0.08 and p = 0.05 (i.e., below the expected 
position of the three-critical fixed point, px—O-l, Fig. Q]) 
and for system sizes up to L — 1024. It is seen that, when 
oscillations corresponding to the lowest discrete states are 
discarded, data for different L nicely combine in a single 
smooth curve corresponding to the thermodynamic-limit 
DoS. At moderately low energy, this curve is well fitted 
Eq. ([3]). So, up to this point, the behavior appears to be 
consistent with the expectation that the RG flow is di- 
rected toward the clean Ising fixed point. However, below 
E~10- 2 {E~10- 3 ) for p=0.08 (p=0.05), the DoS satu- 
rates. For p=0.08, it shows even an upturn for E < 10~ 3 ; 
curves for p—0.05 suggest a similar tendency. 

We do not have an unambiguous interpretation of these 
surprising findings and can only speculate about possible 
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scenarios. 

(i) One possibility is that the position px of the tri- 
critical point T is, in fact, not px — 0.1 as was found in 
Ref. [23J (the phase diagram is reproduced in our Fig. [T]) 
but rather considerably smaller, Pt< 0.05. In addition 
to a conflict with the data of Ref. (23, it would be quite 
surprising if such a numerically small value of px should 
arise. Also, it is pretty unexpected that after having 
reached the value of DoS as low as g ~ 10~ 3 (lower panel 
of Fig. [7]), the system flows toward the metallic fixed 
point. 

(ii) A more sophisticated scenario that would allow rec- 
onciliation of our results with those of Ref. [HI is that, in 
fact, there are two fixed points on the expected quantum 
Hall transition line sin 2 a = 1/2. Namely, in addition to 
the tricritical point px, there is a repulsive fixed point at 
some pn < px- This point would then act as a "flow split- 
ter" which is similar to the role of the Nishimori point in 
the RBIM (hence, the subscript N). 

Taking this idea further, we could imagine that the 
tricritical point T also has a counterpart in the RBIM, 
namely, the zero-temperature transition point from the 
ferromagnetic phase into the spin glass. Furthermore, 
it is possible that the metallic phase of CF model con- 
nects to the spin-glass line of RBIM that exists at zero 
temperature for sufficiently strong disorder. 

(iii) Another nontrivial possibility is that the RG treat- 
ment of the theory of Dirac fermions with Gaussian ran- 
dom mass is, in fact, insufficient, and some effects - pos- 
sibly of nonperturbative origin - eventually drive the sys- 
tem away from the clean Ising fixed point. 

In order to decide which of these scenarios take place, 
further work (analytical as well as numerical) is appar- 
ently needed. We will return to possible directions of 
future research in Sec. HVl 



IV. CONCLUSIONS 

In summary, we have presented a numerical investi- 
gation of the density of states g(E) in the Cho-Fishcr 
network model. The model is a generic two-dimensional 
representative of the symmetry class D describing dis- 
ordered superconductors with broken spin-rotation and 
time-reversal invariances, and shows the thermal quan- 
tum Hall effect. At a sufficiently large concentration p of 
defects, the DoS has a logarithmic divergence as E — > 
with superimposed random-matrix-theory oscillations, in 
agreement with analytical predictions for the thermal- 
metal phase, given in Eqs. ^ and 0. Reducing p, 
we find a transition into localized phases (insulator and 
quantized Hall conductor) with g(E) finite at E=0. 

At the plateau transition between these phases, the 
DoS tends to vanish as E — > in agreement with the be- 
havior g(E) ~ \E\ \n(l/\E\) derived from the theory for 
Dirac fermions with random mass, Eq. ([3]). However, at 
lowest E, this behavior breaks down, and DoS saturates 
and even shows an upturn. We do not have an unambigu- 



ous explanation for this behavior; more work is needed in 
order to understand better the properties of the system at 
the quantum Hall transition line. We imagine an analysis 
of a phase diagram of a wider family of class-D models, 
see Ref. |35j. One can also study the dependence of re- 
sults on the microscopic model in numerical simulations. 
In particular, a possible generalization of the Cho-Fisher 
model is presented in Appendix [Bj Further, one should 
study different observables that are more susceptible to 
the critical behavior and would give additional informa- 
tion about the system at the expected critical line. In 
particular, an analysis of the wave-function statistics is 
expected to be useful in this respect^ 
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APPENDIX A: RENORMALIZATION GROUP 
FOR DIRAC FERMIONS WITH RANDOM MASS 

In this Appendix, we sketch the RG analysis of the 
DoS of Dirac fermions with random mass, leading to 
Eq. © at criticality and Eq. © away from the criti- 
cal line. Our presentation largely follows Ref. [2l| though 
we depart from that work at the end. The Hamiltonian 
has the form 



H 



-id x a x - idyOy + m(r)a z , 



(Al) 



where cr M are Pauli matrices. It satisfies the Hermiticity 
H = W and the class-D symmetry, H = —<r x H T <j x . The 
mass m(r) is a Gaussian random variable, with 



(m(r)m(r')) = 2g M S(r - r'). 



(A2) 



Introducing the field-theoretical representation, perform- 
ing the disorder averaging, and carrying out the RG anal- 
ysis (see Ref. [2l|), one gets in one-loop order the scaling 
equations for renormalization of the disorder strength g^j 
and the energy E, 



dg 



M 



2g 2 



d\nl 
dE 

dha 



M . 



7T / 



(A3) 
(A4) 



Here, i is the ultraviolet cutoff. The model is originally 
defined with a microscopic cutoff £q (lattice spacing), so 
that E = E(£ ) and gu = 9m(4)- Our aim is to analyze 
the DoS, g{E) = g(£(4),jM^o),4)' Integrating these 
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crossing over into Eq. (|A8|) . is 



L 



l + 2 ^ln^ 



-1/2 



(A9) 



Up to now, we considered a system exactly at criti- 
cality (a = 7r/4 in the network model). Moving slightly 
off criticality, we enter the localized phase, with a lo- 
calization length £ < oo. Since the states in different 
localization volumes are essentially independent, the lo- 
calization length will play the same role as the system 
size L in Eqs. and (|A9|) . Substituting £ for L, we 
obtain Eqs. ([8]) and ([9]) of the main text. 



APPENDIX B: 



CONTINUOUS-DISORDER 
MODEL 



FIG. 8: DoS near E — on the self-dual line ao = n/4 for 
the CDM with plaquette coupling parameter a continuously 
distributed within the window Aa. Evolution from metallic to 
critical behavior is clearly observed with decreasing disorder 
strength. The four data curves (from top to bottom) are for 
Aa = 1.0 (x), 0.8 (o), 0.6 (□), and 0.4 (o) with system size 
L=128. Inset: g(E)/\E\ on the log-linear scale (see inset of 

Fig. [Hi. 



equations, one gets 



JmW = 9m(*o) + - In T , 



E{£) = E(£ ) 



2g M {£ ) £ 



1/2 



(A5) 



(A6) 



yielding 



(E(£),g M (£),£) = g(E(£ ),g M (£ )J ) 



1 



■ ml 

£o 



1-1/2 



(A7) 



There are two ways in which the renormalization can 
terminate: (i) E{£) reaches the bandwidth (~ 1) or (ii) 
£ reaches the system size L. In the first case we use the 
fact that for large E(£) the DoS is essentially unaffected 
by disorder, g(E(£)) ~ \E(£)\/2w, yielding Eq. © of the 
main text for \E\ > Eq, where E is given by Eq. (|A9[) . 
In the second case, we get after renormalization a system 
with just few degrees of freedom (ultraviolet cutoff of the 
order of the system size L) and a small energy E(L) <C 1. 
Since the DoS is finite at E=0 in the random-matrix 
theory of class D, we find that 



9(E) 



to 



2^1n^ 



1 1/2 



E <C Eq. 



(A8) 



Comparing Eqs. and (|A9|) . we find that the energy 
-Eo at which the behavior ^ saturates in a finite system, 



The CF model can be defined as a network model with 
the following distribution of plaquette couplings: 



V{a) = (l-p)tf(a- 



V P 

■ a o) + ^ S ( a + a o) + 2 5 ( a + "0 -7TJ- 

(Bl) 

It is natural to expect that properties will be qualitatively 
the same for any model with a generic distribution V(a). 
A precise definition of the word "generic" is far from 
trivial in the present case. In particular, we know that 
the distribution 



V(a) = (1 - p)S(a - ao) + pS(a + a Q ) 



(B2) 



corresponds to RBIM which does not possess a metallic 
phase and thus is not generic. 

In order to test the expectation of (restricted) uni- 
versality, we define a continuous-disorder model (CDM), 
with a Gaussian distribution for the angle a. The cen- 
ter ao of the distribution determines the breaking of the 
symmetry between the C and S plaquettes and thus gov- 
erns the plateau transition. The width Acv determines 
the strength of disorder and therefore replaces the pa- 
rameter p of the Cho-Fisher model. 

In Fig. [51 we show the evolution of DoS on the self-dual 
line, ao = 7r/4, with disorder strength decreasing from 
Aa = 1.0 to 0.4. We see that the behavior is very similar 
to that in the Cho-Fisher model, Fig. [SJ we observe a 
transition from a metal (divergent DoS) to the critical 
region. Furthermore, the families of the curves look es- 
sentially identical in both cases. We interpret this as a 
confirmation of the generic character of the Cho-Fisher 
model. 

We have not studied the behavior of DoS at weak dis- 
order for lowest energies (where we observed a surprising 
upturn in the case of Cho-Fisher model), relegating a 
detailed investigation of this problem to future workj2£ 
where not only spectral properties but also those of wave 
functions will be analyzed. 
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